%% Calculate signal lifetime from trachsFinal data
clear;

[filename pathname] = uigetfile('.mat', 'select result from thunderstorm processing','multiselect','on');
for ii = 1:length(filename)
    load(filename{ii});
    for pp = 1:length(tracksFinal)
        Test = find(tracksFinal(pp).Coord(1,1)==1);
        
        if Test == 1
            Lifetime(pp) = length(tracksFinal(pp).Coord(:,1));
            Intensity{pp} = tracksFinal(pp).Coord(:,5);
        end 
    end
        
        lifetime_spot.lifetime = Lifetime;
        lifetime_spot.Intensity = Intensity;
        filecurr = filename{ii};
        fileroot = filecurr(1:find(filecurr=='.')-1);
        filesave = ['Lifetime-' fileroot '.mat'];
        save([pathname filesave],'lifetime_spot');
    
    clear Spots SpotsLink tracksFinal Lifetime Intensity lifetime_spot
    
end

%% Check intensity histogram
clear;
[filename pathname] = uigetfile('.mat', 'select result from thunderstorm processing','multiselect','on');
for ii = 1:length(filename)
    load(filename{ii});
    for nn=1:length(lifetime_spot.Intensity)
    Int(nn) = lifetime_spot.Intensity{nn}(1,1)
    end
    [pathstr,name,ext] = fileparts(filename{ii});
    Intensity_measurement(ii).filename = [name];
    Intensity_measurement(ii).Intensity = Int;
    clear Int
end
save('Intensity_measure.mat','Intensity_measurement');

%% make a histogram
    x1 = [];
    x2 = [];
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    hty = [1:300:10000];
    x1 = Intensity_222fMp53;
    h1 = histogram(x1,hty,'normalization','probability');
    xlabel('Intensity (photons)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'Hist_222fM_p53.png','-dpng','-r300'); 
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    x2 = No_p53_Nodetector_intensity;
    h2 = histogram(x2,hty,'normalization','probability');
    xlabel('Intensity (photons)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'Hist_No_p53_Nodetector.png','-dpng','-r300'); 
    
    
    mean_Nop53nodetector = mean(No_p53_Nodetector_intensity);
    std_Nop53nodetector = std(No_p53_Nodetector_intensity);
    
    %% Remove autofluorescence spots from p53 sample
    clear;

[filename pathname] = uigetfile('.mat', 'select result from thunderstorm processing','multiselect','on');
for ii = 1:length(filename)
    load(filename{ii});
    for pp = 1:length(tracksFinal)
        Test = find(tracksFinal(pp).Coord(1,1)==1);
        
        if Test == 1
            if tracksFinal(pp).Coord(1,5) > 1090; % mean of autofluorescence spot+std of autofluorescence spot 
            Lifetime(pp) = length(tracksFinal(pp).Coord(:,1));
            Intensity{pp} = tracksFinal(pp).Coord(:,5);
            else
                Lifetime(pp) = NaN;
                Intensity{pp} = [];
            end 
        end 
    end
        
        lifetime_spot.lifetime = Lifetime;
        lifetime_spot.Intensity = Intensity;
        filecurr = filename{ii};
        fileroot = filecurr(1:find(filecurr=='.')-1);
        filesave = ['A_Lifetime-' fileroot '.mat'];
        save([pathname filesave],'lifetime_spot');
    
    clear Spots SpotsLink tracksFinal Lifetime Intensity lifetime_spot
    
end
%% make intensity histogram for samples
clear;
[filename pathname] = uigetfile('.mat', 'select result from thunderstorm processing','multiselect','on');
for ii = 1:length(filename)
    load(filename{ii});
    for nn=1:length(lifetime_spot.Intensity)
        if isempty(lifetime_spot.Intensity{nn}) == 0
           Int(nn) = lifetime_spot.Intensity{nn}(1,1);
        else 
           Int(nn) = NaN;
        end 
    end
    [pathstr,name,ext] = fileparts(filename{ii});
    Intensity_measurement(ii).filename = [name];
    Intensity_measurement(ii).Intensity = Int;
    clear Int lifetime_spot
end

%% Combine the lifetime from each image of the same condition
load('Lifetime_data_No p53- 50%plasmaNo detector.mat');
load('Lifetime_data_222fM p53- 50%plasma.mat');
load('Lifetime_data_adjusted_222fM p53- 50%plasma.mat');
    x1 = [];
    x2 = [];
    x3 = [];
    
    % histogram of autofluorescence signals
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    hty = [50:500:10000];
    x1 = Lifetime_Data_Nodetector*50;
    h1 = histogram(x1,hty,'normalization','probability');
    xlabel('Lifetime of spot (ms)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'lifetime_no_detector.png','-dpng','-r300'); 
    
    % histogram of specific signals with autofluorescence spots
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    x2 = Lifetime_Data_222fMp53*50;
    h2 = histogram(x2,hty,'normalization','probability');
    xlabel('Lifetime of spot (ms)','FontSize',14);
    ylabel('Probability','FontSize',14);
    ylim([0 0.2]);
    print(gcf,'lifetime_222fM_p53.png','-dpng','-r300'); 
    
    % histogram of specific signals without autofluorescence spots
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    x3 = Lifetime_data_adjusted_222fMp53*50;
    h3 = histogram(x3,hty,'normalization','probability');
    xlabel('Lifetime of spot (ms)','FontSize',14);
    ylabel('Probability','FontSize',14);
    ylim([0 0.2]);
    print(gcf,'lifetime_adjusted_222fM_p53.png','-dpng','-r300'); 
    
    lifetime_summary(1,1) = mean(x1);
    lifetime_summary(1,2) = std(x1);
    lifetime_summary(2,1) = mean(x2);
    lifetime_summary(2,2) = std(x2);
    lifetime_summary(3,1) = mean(x3,'omitnan');
    lifetime_summary(3,2) = std(x3,'omitnan');
    
    save('lifetime_summary.mat','lifetime_summary');
    
    
    
    
%% make a histogram of non-specific binding signals and autofluorescence signals
   
    x1 = [];
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    hty = [1:300:10000];
    x1 = Intensity_Nop53;
    h1 = histogram(x1,hty,'normalization','probability');
    xlabel('Intensity (photons)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'Hist_No_p53.png','-dpng','-r300'); 
    
    %% make a histogram of non-specific binding signals 
 
    x1 = [];
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    hty = [1:300:10000];
    x1 = Intensity_Nop53;
    h1 = histogram(x1,hty,'normalization','probability');
    xlabel('Intensity (photons)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'Hist_adjusted_No_p53.png','-dpng','-r300'); 
 %% make a intensity histogram of 75fM p53 sample
    x1 = [];
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    hty = [1:300:10000];
    x1 = Intensity_75fMp53;
    h1 = histogram(x1,hty,'normalization','probability');
    xlabel('Intensity (photons)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'Hist_adjusted_75fM_p53.png','-dpng','-r300'); 
    
    %% make a intensity histogram of 8fM p53 sample
    x1 = [];
    figure;
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    hty = [1:300:10000];
    x1 = Intensity_8fMp53;
    h1 = histogram(x1,hty,'normalization','probability');
    xlabel('Intensity (photons)','FontSize',14);
    ylabel('Probability','FontSize',14);
    print(gcf,'Hist_adjusted_8fM_p53.png','-dpng','-r300'); 
    
%% Separate non-specific binding signals and autofluorescence signals
clear;
[filename pathname] = uigetfile('.mat', 'select result from thunderstorm processing','multiselect','on');
for ii = 1:length(filename)
    load(filename{ii});
    for pp = 1:length(tracksFinal)
        Test = find(tracksFinal(pp).Coord(1,1)==1);
        
        if Test == 1
            if tracksFinal(pp).Coord(1,5) > 1090; % mean of autofluorescence spot+std of autofluorescence spot 
            Lifetime(pp) = length(tracksFinal(pp).Coord(:,1));
            Intensity{pp} = tracksFinal(pp).Coord(:,5);
            else
                Lifetime(pp) = NaN;
                Intensity{pp} = [];
            end 
        end 
    end
        
        lifetime_spot.lifetime = Lifetime;
        lifetime_spot.Intensity = Intensity;
        filecurr = filename{ii};
        fileroot = filecurr(1:find(filecurr=='.')-1);
        filesave = ['A_Lifetime-' fileroot '.mat'];
        save([pathname filesave],'lifetime_spot');
    
    clear Spots SpotsLink tracksFinal Lifetime Intensity lifetime_spot
    
end

%% make a lifetime histogram of non-specific binding signals
 % histogram of specific signals without autofluorescence spots
    x1 = [];
    figure;
    hty = [50:500:10000];
    ax1 = axes('Position',[0.3 0.3 0.3 0.4]);
    ax.ActivePositionProperty = 'outerposition';
    x1 = Lifetime_data_adjusted_74fMp53*50;
    h3 = histogram(x1,hty,'normalization','probability');
    xlabel('Lifetime of spot (ms)','FontSize',14);
    ylabel('Probability','FontSize',14);
    ylim([0 0.2]);
    print(gcf,'lifetime_adjusted_74fM_p53.png','-dpng','-r300'); 
    
    d = mean(x1,'omitnan');
    
    